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It has been shown experimentally that contact interactions may influence the migration of cancer 
cells. Previous works have modelized this thanks to stochastic, discrete models (cellular automata) 
at the cell level. However, for the study of the growth of real-size tumors with several million cells, 
it is best to use a macroscopic model having the form of a partial differential equation (PDE) for 
the density of cells. The difficulty is to predict the effect, at the macroscopic scale, of contact 
interactions that take place at the microscopic scale. To address this, we use a multiscale approach: 
starting from a very simple, yet experimentally validated, microscopic model of migration with 
contact interactions, we derive a macroscopic model. We show that a diffusion equation arises, as 
is often postulated in the field of glioma modeling, but it is nonlinear because of the interactions. 
We give the explicit dependence of diffusivity on the cell density and on a parameter governing 
cell-cell interactions. We discuss in detail the conditions of validity of the approximations used in 
the derivation, and we compare analytic results from our PDE to numerical simulations and to some 
tn vitro experiments. We notice that the family of microscopic models we started from includes as 
special cases some kinetically constrained models that were introduced for the study of the physics 
of glasses, supercooled liquids, and jamming systems. 

PACS numbers: 87.18.Gh, 87.10.Ed, 05.10.-a, 87.10.Hk, 87.19.xj, 87.19.1k 



I. INTRODUCTION 

The migration of tumor cells plays a principal role in 
tumor malignancy, particularly for brain tumors. The 
fact that glioma cells migrate fast contributes greatly 
to glioblastoma lethality and thwarts therapy strategies 
based on tumor resection Understanding the mech- 
anisms of migration can be of utmost importance when 
it comes to devising efficient treatments 0, @| • In^ vitro 
studies combined with mathematical models constitute a 
promising exploration path, a fact supported by several 
high-quality results [1, Q . 

How does one go about modeling cell migration? There 
are essentially two approaches: on one hand, "macro- 
scopic" models that take the mathematical form of one 
or several partial differential equations (PDE) 0, Q, 
and, on the other hand, "microscopic" models where in- 
dividual cells are represented and evolve according to 
some (stochastic) rules: the so-called individual based 
models, agent based models, cellular automata, and so 

on 

In a macroscopic model, one eschews all reference 
to the elementary constituent, the cell, and introduces 
macroscopic quantities, like the average cell density, 
which evolve according to a PDE. One expects stochastic 
deviations from the average to be negligible at the scale 
of a large population of cells, provided this population is 
homogeneous (quite often, cell populations in cancer are 
made of cells having distinct pheno types or genotypes, 
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because e.g., of genetic mutations [H, [T^, and this 
should be taken into account in a more general model of 
cancer development). As a matter of fact, the first brain 
tumor models (and the vast majority of models exist- 
ing today, which are essentially refinements of these first 
ones) consider the dynamics of a small number of macro- 
scopic quantities The tumor cell density is, of 
course, incontrovertible, the remaining quantities being 
usually concentrations of chemicals, which, one surmises, 
play a role in the migration. 

Using a PDE has many advantages over a "micro- 
scopic" model. First, even with the ultrasimplification 
of a cellular automaton, one can treat a relatively small 
number of cells, at best a few thousand, a number very far 
from the millions one encounters in a tumor lump. Then, 
there are ready-to-use PDE solvers. Since the PDEs are 
deterministic, one simulation is enough where one should 
perform many simulations of the cellular automaton to 
average over the stochastic noise. The number of dis- 
cretization steps of the space is fixed only by the meso- 
scopic geometry of the medium in which the cells diffuse; 
there is no need to have as many steps as cells. Finally, 
it is easier to get analytical results from a PDE (even 
approximate) than from a cellular automaton. 

In this spirit, the simplest approach to cell migration is 
to assimilate the process to linear diffusion. Then a dif- 
fusion equation may be used to find how the cell density 
evolves with time. 



dp 
dt 



(1) 
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This was, for instance, the approach of Burgess et al. [18| . 
who modeled the evolution of glioblastomas with the help 
of a diffusion equation (plus an exponential growth rep- 
resenting cell proliferation). Refinements of the diffu- 
sion model can be, and have been, introduced. One can, 
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for instance, consider a space-dependent diffusion coeffi- 
cient that varies depending on the nature of the tissue 
in which migration takes place [l^, [IS, HH, [13] ■ Terms 
can be added to Eq. ([T]) representing various chemo- 
tactic effects, either as arbitrary external source terms 
or as quantities that couple Eq. ([T]) to other evolution 
equations [1^ [IJ . When the cells have a velocity or a 
persistence of the direction of motion, one has to study 
(e.g., at short time scales of the cell migration) a different 
type of equatioii, like a hyperb olic PDE or a degenerate 
parabohc PDE ^ [H [lE [H . 

However, modeling of cell migration as sirnple, heatlike 
diffusion (noninteracting random walks [l^HO]) neglects 
important phenomena, even in the absence of chemo- 
taxis: cells cannot penetrate each other, and they often 
interact via contacts with other cells and/or with the 
extracellular matrix. 

But what PDE should one use to take contact inter- 
actions into account in the migration process? It can be 
easier to devise realistic evolution and interaction rules in 
a cell-based model [1, [HI, [3l| . Even though there are situ- 
ations in which it is not possible to discriminate between 
several macroscopic models, especially when only data for 
the evolution of the cell density are at hand [U, [H, [s^l , 
it is of little interest to use a PDE that can later prove 
to be wrong when more data become available and, in 
particular, when one is able to compare predictions for 
the trajectories of individual cells [sl. [Sll. [s^ . 

Sometimes, a microscopic model is derived from a PDE 
through a discretization procedure [34j . But there are 
several ways to discretize space that lead to very dif- 
ferent behavior of the microscopic agents, some of which 
can be unrealistic. (For example, consider the discretized 
population pressure model of [33|: one can show that, 
for some discretization procedure, cells move toward the 
high-density region they are supposed to flee, but for 
some other procedure they move away from it — and 
both correspond to the very same PDE.) At least, one 
needs a biological criterion to choose the "right" dis- 
cretized version of the PDE. 

Therefore, to establish a macroscopic model for tumor 
cell migration with contact interactions via the so-called 
gap junctions between cells, we start from the micro- 
scopic model introduced and successfully compared to 
in vitro migration experiments in [s^. In this setting, 
cell proliferation, apoptosis (cell death), and mutations 
(changes of cell behavior) are negligible (hence the popu- 
lation is homogeneous). The number of cells is constant 
and cells simply move on a collagen substrate. The model 
is defined in Sec. II. (There we also discuss special cases 
in which our model is a kinetically constrained model, 
borrowing some knowledge from the theory of glasses, 
supercooled liquids, and jamming systems.) 

The derivation of the macroscopic model, the so-called 
hydrodynamic limit, is carefully explained in Sec. III. 
The derivation of hydrodynamic limits is a well-known 
technique in mathematical physics [H, 113, [H, [s^ [iOl . 
and it has also been used in mathematical biology [g. 
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field of cell movement and cancer modeling. Since our 
derivation relies on a mean-field like approximation, we 
review in Sec. Ill the circumstances under which the ap- 
proximation fails, and we provide thorough comparisons 
to simulation results. 

The PDE we obtain is a nonlinear diffusion equation. 
Contrary to what is customary in the phenomenological 
models to which we alluded above, the terms of the diffu- 
sion equation depend crucially on the cell-cell interaction 
present in the microscopic model. As our analysis shows, 
these interactions introduce nonlinearities in the diffu- 
sion equation. These nonlinearities are essential in order 
to reproduce the density profiles obtained in in vitro mi- 
gration experiments (Sec. IV). 



II. THE STOCHASTIC CELLULAR 
AUTOMATON MODEL 

In [3^, a cellular automaton (lattice gas) was intro- 
duced to mimic the migration of cancer cells with con- 
tact interactions. The dynamics of migratory motion are 
described through the evolution of points (representing 
the cells) on a grid (or lattice), while the interactions are 
introduced through evolution rules. The impossibility of 
cells to penetrate each other is trivially taken into ac- 
count (no two cells can share the same position on the 
grid), and the scale length of the lattice step is the typical 
size of one cell. 

It is easy to extend such a cellular automaton to cells 
that occupy several sites on the lattice, using e.g., a cellu- 
lar Potts model [lll,lill[il,li3|. But then we would need 
to introduce many parameters (surface tension, binding 
energies, ...), whereas experimental data at hand were 
not sufficient to measure so many values and avoid dubi- 
ous fits. Thus the model of ^35j is an elementary cellu- 
lar automaton with nonextended cells, involving a single 
parameter that takes contact interactions into account. 
This approach has made it possible to model with suc- 
cess the migration of glioma cells over substrates of col- 
lagen or of astrocytes and draw conclusions on the exis- 
tence of homotype and heterotype interactions between 
the cells ^EM- 

In order to define a cellular automaton model, one 
needs two things. First, one must choose the geome- 
try of the lattice. One may work in one, two, or three 
dimensions and consider elementary cells of any shape 
and size. Second, one must specify the update rules of 
the automaton. With the adequate choice of these rules, 
one can mimic the interactions between the elementary 
entities (glioma cells in [s^) which are modeled by the 
automaton. In this paper, we shall follow the prescrip- 
tions we introduced in [35|. We recall them briefly below. 
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A. The geometry 

The cellular automaton introduced in 35] for the de- 
scription of the glioma cell migration is based upon a 
regular hexagonal tiling (or, equivalently, a triangular 
lattice). This choice is dictated by the fact that the 
hexagonal tiling is the less anisotropic among all regu- 
lar tilings of the plane. Thus one can have a very simple 
algorithmic definition of the geometry while keeping the 
constraints due to the symmetry of the lattice at a mini- 
mum. Moreover, in such a lattice each cell is surrounded 
by six sites, a number sufficiently high to allow a cer- 
tain freedom of motion at each update. A generalization 
to three space dimensions is obtained by considering the 
face-centered-cubic (fee) lattice or the hexagonal com- 
pact (he) lattice, where each site has 12 nearest)-neighbor 
sites. It is known that some difhculties arise in three- 
dimensional lattices as compared to two-dimensional lat- 
tices when one aims at finding a cellular automaton, the 
continuous limit of which would be the Navier-Stokes 
equation [i^, but here we have no such constraint: 
instead of starting from a macroscopic equation (such as 
the Navier-Stokes equation) and guessing what cellular 
automaton would mimick it, we start from a given cellu- 
lar automaton and ask for (an approximation to) the cor- 
responding macroscopic equation. The only constraint is 
that the macroscopic equation be a reasonable model of 
the behavior of real cells in experiments. Additionally, 
it turns out that both the fee and he lattices, although 
they do not have the same symmetry properties ^79.] , lead 
to the same macroscopic equation. Therefore, we expect 
that refinements in the choice of the three-dimensional 
lattice would not significantly change our results. 

In each case, we denote by a the distance between the 
centers of two nearest-neighbor sites (lattice step). 

B. The stochastic evolution rules 

Each hexagon can be occupied by at most one cell at 
a time. As a consequence, a cell can move only to a 
free site. The total number of cells is locally conserved 
(but cells may be introduced into or removed from the 
system on the boundaries). Thus our automaton is a 
kind of lattice gas or "box and ball" system [4^ . Notice 
that our model is not a so-called lattice-gas cellular au- 
tomaton (LGCA), where individual particles have both 
a position and a speed: here the cells lose the memory 
of their speed and previous position after each move. In 
the classification of [301 ^ it is a "position jump process," 
not a "velocity jump process." 

At each update of the automaton (discrete time step), 
the process of picking a lattice site at random and apply- 
ing the stochastic evolution rules to this site is repeated 
a number of times equal to the total number of sites, so 
that, on average, each site is updated once at each time 
step [iOl- If the site is empty, nothing happens. If it is 
occupied by one cell, the following rule is applied. It is 



inspired by what we believe is the biological mechanism 
of contact interaction. The interaction paremeter p is a 
fixed, constant number between and 1. We pick at ran- 
dom one of the six neighboring sites, the target site. If 
the target site is occupied, nothing happens. Otherwise, 
the move from the cell to the target site (leaving its de- 
parture site empty) is done with probability p if at least 
one of the two common neighboring sites of the departure 
site and of the target site is occupied, and with probabil- 
ity 1 — p if the two common neighboring sites are empty. 
With these assumptions, p — 1 means that a cell can only 
move if it is in close contact with at least one other cell, 
and it can only move to a position where it will stay in 
contact with a least one of its former neighbors, while, 
for p = 0, a cell moves only in such a way that all con- 
tacts to its former neighbors are broken. If p = i, cells 
are indifferent to their environment (apart from mutual 
exclusion), and we expect that the cell density obeys a 
diffusion equation. Values of p > ^ model situations in 
which cells are reluctant to break established junctions 
with their neighbors, whereas p < ^ involves a kind of 
short-range repulsion but only between cells that were 
in close contact, unlike what is usually called repulsion 
between, e.g., electric charges. If these rules seem odd 
from a physical point of view (they cannot be related 
to some energy or free energy and convey a strange no- 
tion of repulsion, and in particular they do not reduce 
to a special case of the model of 0]), they incorporate 
the biologically reasonable assumption that a cell makes 
its decision regarding where it goes by considering only 
information it can get from its immediate vicinity, e.g., 
through established communication junctions with the 
cells with which it is in contact. 

The preceding rules can readily be extended to the 
three-dimensional hexagonal compact or fee lattices; the 
only change is that a pair of a departure and a target site 
has four neighbors in common instead of two. 

The setting considered in [35|] consists in a central part 
of the lattice occupied by the equivalent of the glioma cell 
spheroid, which, we assumed, acts as a cell source. Once 
a free position in the hexagons surrounding the center 
appears, it is immediately occupied by a cell "released 
from the center." In [ssj . we have explained how one can 
(and must) calibrate the model so as to allow quantitative 
comparisons with experiment. 



C. Kinetically constrained models 

Our cellular automaton contains as special cases (for 
p = 1 and 0) two kinetically constrained models (KCMs). 
These models were introduced in statistical physics to 
mimic the dynamics of glasses, supercooled liquids, and 
jamming systems j49l[50|. 

In a KCM, a microscopic object (in our context, cells) 
can move to a free position only if some requirement 
about its neighborhood is fuUfiUed (e.g., only if it has 
at least another empty position around itself). As a con- 
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sequence, at least for some concentrations of surrounding 
cells, the motion of an individual cell is not a Brownian 
motion, as it is in a diffusive model such as the simple 
symmetric exclusion processes (p close to i). However, 
if one does a "zoom out" or coarse-graining of the trajec- 
tory of an individual cell, using characteristic length and 
time scales (which depend on the concentration of cells) 
as new unit length and time, the result looks Brown- 
ian [5l|. Above these scales, a diffusion equation should 
hold. 

The origin of this separation of scales and the existence 
of the characteristic length and time scales is easy to un- 
derstand for the so-called noncooperative KCMs (mod- 
eling strong glass- former materials). There, cells are 
blocked unless they are close to some mo del- dependent 
defect or excitation that has a finite size. These exci- 
tations travel across the system in a kind of diffusive 
motion, and the configuration of the cells in a region will 
be changed or relaxed only after the region has been tra- 
versed by an excitation. Even though the excitations, 
and the individual cells while they are close to an exci- 
tation, diffuse "quickly," a cell goes "slowly" away from 
its initial position because its trajectory is an alterna- 
tion of long nonmoving phases (when the cell is away 
from any excitation) and short diffusive phases (when 
the cell is close to an excitation) [5l|, [s^l • The character- 
istic time scale is the typical time between two stays of 
an excitation in the neighborhood of a cell, and the char- 
acteristic length is the typical size of the region explored 
by a cell during such a stay. Over times much greater 
than the characteristic time, the trajectory of a cell is 
made of many nearly independent non-moving/quickly 
diffusive patterns and, according to the central limit the- 
orem, the cell undergoes a Brownian motion with an ef- 
fective diffusion constant. Pick's law is obeyed and the 
coarse-grained model is like a lattice gas without kinetic 
constraints [soj . 

In the case of the cooperative KCMs (modeling fragile 
glass-formers) , no traveling finite-size excitation facilitat- 
ing the motion of cells exists. Instead, the moves of the 
cells are organized hierarchically and form structures that 
can be infinitely large [s^]: some changes of the system 
require the participation of an extensive fraction of the 
cells, hence the term "cooperative". A complete block- 
ing at a finite concentration of cells can even be observed 
for some models [53| . Still, in the absence of such block- 
ing, it is possible to define a characteristic time and a 
characteristic scale above which a diffusive behavior is 
recovered and Pick's law is obeyed [54j . 

In our cellular automaton, both the noncooperative 
and cooperative behaviors take place, for p = 1 and 0, 
respectively (and the behavior for p close to 0, on the one 
hand, and for p close to 1, on the other, is qualitatively 
the same as in these two limit situations). 

Por p = 1, our model is a noncooperative KCM: ac- 
cording to the rules of the automaton, an isolated cell 
cannot move, but a cluster of two neighboring cells has 
a finite probability to go anywhere in the lattice and it 



can be "used" to move any isolated cell (it plays the role 
of the aforementioned excitations) . 

Por p = 0, our model is a cooperative KCM: this is 
the very same model as the "2-triangular lattice gas" or 
(2)-TLG studied in the context of the slow dynamics of 
glasses [53,[5i,[5l[5i|. 

To conclude, for general p, our model interpolates be- 
tween a cooperative KCM (p = 0), a simple symmetric 
exclusion process (p = 5) and a noncooperative KCM 
(P=l)- 



III. THE HYDRODYNAMIC LIMIT 

In this section, we explain in detail how we derive 
an approximate macroscopic, deterministic model (which 
takes the form of a partial differential equation) from the 
microscopic, stochastic cellular automaton. This kind of 
technique has been used in physics since the 1980s to take 
the so-called hydrodynamic limit of a stochastic model. 
We recall the basics for readers who are not familiar with 
it. Other readers may be interested in our discussion of 
rigorous results and of the quality and limits of our ap- 
proximation. 

The hydrodynamic limit exists for some model defined 
on a lattice if, as the size of the lattice goes to infinity, 
one can cut the lattice into parts such that (i) the parts 
are negligibly small with respect to the whole lattice, 
(ii) the values that macroscopic quantities (like the cell 
concentration) take for each part have vanishingly small 
fluctuations, and (iii) the averages of these values obey 
some partial differential equation(s). This process is a 
kind a coarse-graining. Similarly, water, though made of 
discrete molecules, looks like a fluid composed of homo- 
geneous droplets, because each droplet contains a huge 
number of molecules that are typically spread homoge- 
neously in the droplet rather than all packed in one-half 
of the droplet. 

The establishment of a PDE for macroscopic quan- 
tities after coarse-graining has been extensively studied 
and can be made mathematically rigorous in a number 
of cases (see, e.g., [H, S HI, iSl)- We would like to 
show that, although one might think that this continu- 
ous time and continuous space approximation is useful 
only in the case of very large lattices, it yields results 
in excellent agreement with the results from our cellular 
automaton for lattices as small as 16x16. Indeed, even if 
in one simulation (or in one experiment) the microscopic 
concentration of tumor cells is discrete (it can take only 
two values on each lattice site, depending on whether the 
site is full or empty), the average of this concentration 
over independent simulations of the cellular automaton 
(or over independent experiments) is well predicted by 
the PDE. Purthermore, there are many ways to perform 
the approximation, some of which being rather involved, 
but we show that, here, an elementary procedure that can 
easily be applied to other cellular automata is sufficient. 



5 



A. On the existence and scale of the validity of the 
hydrodynamic limit 

Several authors have proven rigorously that some cellu- 
lar automata have a nonlinear diffusion equation as their 
macroscopic scaling limit [13, [H, HI, In particular, 
the theorem of fEd] proves that our model has a hydrody- 
namic limit for all values of p between and 1 excluded 
(we have to exclude and 1 because, for them, some 
transition rates are and the theorem does not apply). 
However, it seems rather complicated to get an explicit 
formula for the corresponding nonlinear diffusion equa- 
tion from this theorem. On the contrary, our computa- 
tion provides approximate but explicit formulas that are 
very useful most of the time. 

When p = 1 or 0, one needs to be more careful: then, 
we have seen in Sec. II that the cellular automaton is a 
KCM where Pick's law is obeyed only above some char- 
acteristic length and time scales. 

Very recently, Gongalves et al. constructed some 
KCM, the hydrodynamic limit of which would be the 
porous media equations, and they were able to prove this 
rigorously [6l| . For p = 1, our model is only slightly dif- 
ferent from theirs and we believe that generalizing their 
proof to ours would be possible. This is an indication 
that (nonlinear) diffusion should hold above some spa- 
tial scale, but in practice it is necessary to know how 
large this scale is. We found that, for p = 1, and more 
generally for all values of p not close to zero, this scale 
is of the order of the elementary lattice spacing — the 
agreement between simulations of the cellular automaton 
and solutions of the PDE is excellent already at this scale 
(see later). 

For p — 0, we know from the previous section that our 
model is a cooperative KCM (the "2-triangular lattice 
gas" i56i]). To the best of our knowledge, there is 
no rigorous proof of the hydrodynamic limit in this case 
or in a similar case, but it is highly probable that this 
limit exists because the length and time scales we have 
discussed above exist in particular for this model [131 
(they diverge as the cell concentration goes to 1, but 
are otherwise finite). We found that our approximation 
breaks down as p gets close to 0, which is a sign that the 
kinetics cannot be understood only by looking at one cell 
and its nearest neighbors. Still, p needs to be as small as 
0.05 for one to see a significant discrepancy between our 
approximation and the automaton. 



B. Formalism and computation 

We now proceed with our computation. On the basis of 
the previous considerations, this approximate computa- 
tion should be valid only if the spatial scale above which a 
macroscopic, coarse-grained diffusive behaviour is of the 
order of the lattice size — fortunately, this is the case for 
most values of p. 

Our procedure involves first taking averages over the 



stochastic noise, neglecting correlations, in order to yield 
deterministic equations on the lattice, then taking the 
continuous space limit to go from the microscopic to the 
macroscopic scale. Of course, averaging over the noise 
amounts to losing information about the true process. 
Describing the behavior of a few cells on the lattice by 
the sole average of the cell concentration on each lattice 
site would be a crude approximation, of little interest. 
However, one can hope that this approximation is not so 
bad in the presence of a large number of cells, and that 
the typical fluctuations of interesting macroscopic quan- 
tities such as the cell concentration become very small as 
the system grows. Conversely, large fluctuations of the 
macroscopic quantities would become very rare as the 
system grows. As we shall see, this does happen for our 
model and the deterministic, continuous space aproxima- 
tion is already excellent for rather small systems. 



1. Averaging over the stochastic noise 

To perform the averaging and neglect correlations, 
we use a technique that is similar to the Hubbard- 
Stratonovich transform and which was already used 
to study another cellular automaton, the contact pro- 
cess [g^l- This technique may be used to build a pertur- 
bative expansion. Here, we will actually stop at the first 
order, since it yields results that are already in excel- 
lent agreement with our numerical simulations. There- 
fore, introducing the technique is not necessary and one 
can get the same results using other means. But we 
explain it because its formalism is convenient and com- 
pact, especially if one wants to use computer algebra to 
compute the macroscopic PDE on complicated or high- 
dimensional lattices. In addition, it yields access to sys- 
tematic perturbative corrections (we can, in principle, 
get closer and closer to the exact result by keeping more 
and more terms) on one hand, and to large deviations 
properties, on the other hand (see [6^ for an example 
where large deviations properties are crucial). 

Using the now standard "quantumlike" formalism for 
master equations, introduced a long time ago [63, 64, 6^, 
l66j , we start by representing the probability distributions 
of the configurations of the cellular automaton as vectors 
in a simplex. The stochastic rules are encoded as a linear 
operator acting on this space. 

For each site i of the lattice, we introduce two basis 
vectors jOj) and to encode the situations in which 
the site is empty and occupied by one cell, respectively. 
If the lattice had only one site, say i, only two config- 
urations would be possible ("the site is occupied by a 
cell" or "the site is empty") and it would be enough 
to specify the occupation probability p of the (unique) 
site to characterize entirely the statistical distribution of 
the configurations under some fixed environment. In our 
formalism, we represent this distribution by the vector 
Vi := p -I- (1 — p) \0i); the coefficient of each basis vec- 
tor is the probability of observing the corresponding con- 
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figuration (there is no need to take the square modulus as 
in quantum mechanics). If the lattice is composed of two 
sites, say i and j, there are four configurations and the 
most general probability distribution can be represented 
by the vector a\0i,0j) + /? |0„ 1^) + 7 Oj) + 6\U,lj) 
with a + (3 + ^ + S = 1 (probabilities must sum up to 
1). Here, the four basis vectors are the tensor prod- 
ucts of the basis vectors for the site i and of the site 
j: \li,Oj) = ® \0j) and so on. In the special case 
in which the sites i and j are statistically independent 
(in which case the probability that i is occupied and j 
is occupied is the product of the occupation probabilities 
of i and j, say pi and pj, respectively), this vector fac- 
torizes as [(1 - pi)\0^) + pi\U)] ® [(1 - pj)\Oj) + Pj\lj)]. 
More generally, if there are N sites, there are 2^ possible 
configurations and our vector space has dimension 2^. 

To explain how we encode the stochastic rules, let us 
take again the simple case in which the lattice has only 
one site. Let us address first the stochastic process of 
"radioactive decay": when the site is full, it has a con- 
stant probability per unit time a to become an empty 
site, and, when it is empty, it remains so. The occu- 
pation probability pi{t) of the site obeys the differential 
equation 

dpjdt = -ap,{t). (2) 

This equation is called the master equation or the 
Chapman-Kolmogorov equation. The vector that en- 
codes the probability distribution, Vi{t), obeys 

dvi{t)/dt = a(ci - cfci)vi(t), (3) 

where q and are 2x2 matrices: in the basis |0), |1), 




The symbols for these matrices are borrowed from quan- 
tum mechanics, since c and behave like annihilation 
and creation matrices, respectively (for so-called hard- 
core bosons). The evolution operator Ci — cfci of vector 
Vi has two terms: the first one, Ci, lets the probability 
of the configuration |0) increase proportionnally to the 
probability of the configuration |1), and, without the sec- 
ond one, probability conservation would not be ensured 
because the probability of the configuration |1) would be 
unchanged. Let us give a second example in which two 
sites i and j are coupled: consider the stochastic process 
in which a particle is transferred from the site i (if site 
i is occupied) to the site j (if it is free) at rate k, i be- 
coming empty and j full. The evolution equation of the 
vector V encoding the probability distribution of the 4 
configurations of the couple of sites reads 

dv/dt = k[c+q - (I - c+Cj)c+c,]v{t). (5) 

Here we extend the definition of the matrices. For in- 
stance, we understand as a 4 x 4 matrix, acting on the 
whole space of dimension four where v lives, equal to the 



tensor product of the former matrix Ci (which acts on 
the two-dimensional space of Vi) and of the identity ma- 
trix of the two-dimensional space of Vj . I is the identity 
matrix acting on the whole space. The first term of the 
evolution operator, c^Q, lets the probability of occupa- 
tion of j increase proportionally to the probability that i 
was full and that j was empty. The second term ensures 
overall probability conservation. 

It is easy to generalize this to any stochastic process. If 
the process has several rules that apply in parallel (e.g., 
"for any couple of neighboring sites i and j, a particle 
can jump from i to j if i is full and j empty, making i 
empty and j full"), the evolution matrix is the sum of 
the evolution matrices encoding each individual rule [in 
the example above, it would be the sum over the couples 
of neighboring sites i and j of the matrices c~''Ci — (1 — 
c1j Cj)cl Ci\. It is equivalent to defining a process by its 
stochastic rules and by the expression of the evolution 
matrix of its configuration probability vector v. 

Applying this technique, we find that our cell migra- 
tion process has the following evolution matrix: 

W=- ^ [c+Cj-(l-C+Cj)c+Ci] F ({c^Cfe}fe n.n. of i and j) ■ 
i,j n.n. 

(6) 

"n.n." stands for nearest neighbors on the lattice, and 
z is their number for a given site (coordination number, 
equal to 6 on the hexagonal tiling and to 12 on the three- 
dimensional fee or he lattices). F is a polynom of two or 
several variables; for each i and j, it is applied to the 
matrices c^Ck, cf ci and so on where fc, /, ... are the com- 
mon neighbors of i and j on the lattice. The expression 
of F is chosen so that, for each possible configuration of 
the occupation numbers nk, ni, ... of the sites k, I, 
F{nk,ni, . . .) is equal to the transition rate of one cell 
from site i to site j (provided that site i is full and site 
j is empty). It would be easy to study other rules, but, 
to conform to the model defined in the previous section, 
we took in two dimensions on the hexagonal tiling 

F{nk,ni) = p (nk+ni) (3-nA;-n;)/2+(l-p)(l-nfc)(l-n;) 

(7) 

which is equal to 

F{nk,ni) ^ p {nk+ni-nkni) + {l-p){l-nk){l-ni) (8) 

since ~ Uk and nf ~ ni and, in three dimensions on 
the fee or he lattices, 

F{nk,ni,nm,nn) = ps{5 - s)(s^ - 5s + 10)/24-|- 

(1 -p)il- nk){l - ni){l - n„)(l - n„), (9) 

where s is short for Uk + ni + Um + rin. In both dimen- 
sions, we want F to be equal to p if at least one of the 
nearest neighbors of i and j is full, and to 1 — p if all are 
empty. The previous expressions are the simplest ones 
that achieve this. 

We use the technique of Ref. [63] to derive from the 
expression of the evolution operator (W) the mean-field 
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(correlations are neglected) evolution equations for the 
occupation probability Pj it) of each site i. There are 
equivalent techniques (sj, [ill, \^ to get the bare mean- 
field evolution equations of pi{t)] this one can addition- 
ally yield results for the large deviations and systematic 
perturbative corrections. The evolution equation of pi{t) 
in the typical (i.e., most probable) configurations reads 



dtPtit) = -5^^(f)(W^), 



(10) 



where dt denotes time derivative, 9.^.(f) is the functional 
derivative with respect to "^i at time taken at — 0, 
and {W) is the average of the evolution operator. This 
average is obtained by replacing each Ij with 1, each 
cfci with pi{t), each Ci{t) with /9i(i) exp['(/'i(t)], and con- 
sistently each cf{t) with [l — pi{t)] exp[— -0i(t)]. Using the 
expression of W, Eq. ([6]), one gets the evolution equation 
of the typical site occupation probability on the discrete 
lattice. 



dtPt{t) = ^ ^ [pj{t)-pi[t)]F [{pk{t)}k n.n. of i and j] 

(11) 



j n.n. of i 



2. Continuous space limit 

So far, we have obtained N evolution equations for the 
N average occupation probabilities pi [t) of the sites from 
the initial 2^ evolution equations for the probabilities of 
the 2^ configurations of the lattice. We shall replace 
the N quantities pi{t) by a single field p{r,t) such that 
p{ri, t) = pi{t) at all times fi being the position of the 
center of site number i. If the lattice step a is finite, there 
are infinitely many such fields, but to the a — > limit 
there is only one that is regular (twice differentiable). In 
reality, a is of the order of a single cell's diameter and 
it is not vanishing (this is a natural length scale in the 
problem, and it remains relevant after our procedure): 
taking formally the a — > limit means actually that one 
is interested in phenomena that take place on a length 
scale of many individual cells, as is standard in statisti- 
cal mechanics, by restricting to spatial variations of the 
average cell density p that are slow on the length scale a. 

For a given site i, we replace each quantity pj{t) — 
p(fj , t) with the Taylor expansion of p{fj , t) in space 
around r^: introducing the unit vector from i to j, 
:= [fj ~ ri)/a^ we can write 

p{rj,t) = p{ri,t) + aui^j ■Vp{ri,t) + 

{ay 2) u,^j ■ H{n,t)u,^j + 0{a^) (12) 

where the center dot denotes the usual scalar product, 
V is the gradient operator, and H is the Hessian matrix 
of p (here taken at the position fi and at time t). After 
this substitution, all terms of order in a vanish because 
there is locally conservation of the number of cells (no 



proliferation nor apoptosis), and all terms of order a van- 
ish because the evolution rules are symmetric (induce no 
bias between left and right) and the lattice is reflection- 
invariant. Also, the lattice is translation-invariant and 
the sites are coupled all in the same way (to their near- 
est neighbors), thus each evolution equation yields the 
same relation between the field p and its space and time 
derivatives. Since we are interested in the continuous 
space (a — > 0) limit, we disregard terms of order 3 and 
above in a, and after some algebra we find the evolution 
equation of the site occupation probability. 



dtpif, t) = o'dwiD (p(f, t)) Vp(r, t% (13) 



where 



D{p) = (1 - p)/4 + {2p - l)p{l - p/2)/2 



(14) 



on the hexagonal tiling in two dimensions (2D) and 

D{p) (1 - p)/6 + {2p - l)p{A ~Qp + 4p2 - p3)/6 (15) 

on the fee lattice in 3D. For the he lattice, we find the 
very same expression of D{p) as for the fee lattice. The 
expression of D{p) is the value of the function F [Eqs. ([S]) 
and ([9])] when all its arguments are equal to p, up to 
a geometrical factor. In the a ^ limit, the right- 
hand side of Eq. (jl3p vanishes, leaving the trivial equa- 
tion dtp{f,t) ~ 0: this is a sign that nothing interesting 
happens, in the limit of continuous space, at the time 
scale of the updates of individual sites, t. The only time 
scale leading to a nontrivial equation has a unit time of 
the order of 1 /a^ [Slj . This scale is characteristic of diffu- 
sion (in the case of an ad vective phenomenon, we would 
have chosen the time scale 1/a [3^). Equation (fT5|) is a 
diffusion equation with a concentration-dependent diffu- 
sion coefficient D{p). This equation belongs to the family 
of the porous media equations [67[ . A striking difference 
with the heat equation (or with Fick diffusion) comes 
from the possibility that D vanishes at some site occu- 
pation probability p, which may create abrupt fronts in 
the diffusion profile instead of large, Gaussian tails. 

To completely free ourselves from the lattice, which in 
most practical applications does not exist, but was con- 
venient to define the model, we write down the equation 
for the spatial concentration of cells c instead of the site 
occupation probability p. This is done be dividing p by 
the volume of a single site, ^/3a^/2 for the hexagonal 
tiling, and a^/{3V6) for the fee lattice. Because p and 
c are proportional, we can discuss either one in the rest 
of this Section; when we discuss experimental results, we 
shall use only c. 

Before proceeding to a comparison of this result to 
numerical simulations, let us discuss the vector "density 
of cell current," j. We expect it to satisfy Fick's law, with 
a nonlinear diffusivity. To compute J, we first compute 
the net cell current along the lattice link i — > j, i.e., the 
current from site i to site j, 



-4c,) = -F {{p{f,t), . . .})u,^j ■Vpif,t) + 0{a^), 

(16) 
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where r := (r^ + rj)/2 and all arguments of F are equal 
to p(r, t). Then we add the contributions of all links (on 
the hexagonal tiling, there are three types of contribu- 
tions since the links can have three directions), and we 
go to the density of current by multiplying each contri- 
bution by the density of links with that direction on a 
hyperplane orthogonal to that direction. On the hexag- 
onal tiling, if some links are, say, vertical, then we count 
the density of intersections of a horizontal line with links. 
These intersections have a periodic pattern and can be 
grouped four by four; a group has the length V3a and 
contains two vertical links (each of these contributes one), 
and two links making an angle 7r/3 with the vertical di- 
rection [each of these contributes cos(7r/3)], hence the 
factor \/3/a- We find on the hexagonal tiling 

jlr, t) = -^DW, i)]Vp(r, t). (17) 

Putting this expression into the conservation equation 
for the number of cells (there is neither apoptosis nor 
proliferation in our model), 

9tc(f,i) = -divj(f,i), (18) 

and using that c = 2(0/(\/3a^) (on the hexagonal tiling), 
we find again Eq. (fT3|) . 

C. Comparison to steady-state simulations of the 
automaton 

In order to test the analytical results p3|) and p4)) . 
we did some simulation of the cellular automaton in a 
simple geometry. More realistic simulations are the topic 
of the next Section. Since, as will be discussed later on, 
the equilibrium state of the cellular automaton is trivial 
(all configurations of the tumor cells are equally proba- 
ble), we chose a setting where a nonequilibrium steady 
state of the tumor cells can exist: a cylinder made of 
riL + 1 rings of cells (its total length is (n^ + l)a\/3/2) 
and base circumference H = nua connected to a reser- 
voir full of cells {p — 1) at the left end and to an empty 
reservoir {p — 0) at the right end, as shown in Fig. [T] 
{riL and nu are two integers). The first and the last 
ring of cells of the cylinder belong to the full and to the 
empty reservoirs respectively; hence the effective length 
over which diffusion takes place is L := nLa\/Zl2. We 
chose a cylinder rather than a simple rectangle in order 
to make the boundary effects in the direction perpendic- 
ular to the flow small, if not vanishing. We introduce 
two coordinates on the cylinder: x is measured in the 
direction parallel to the axis and y along the shortest cir- 
cles. To simulate the reservoirs, on the boundaries of the 
system, cells behave under special evolution rules: each 
time one of the sites on the left border of the lattice is 
left by a cell, it becomes occupied again by a new cell 
with no delay, and each time a cell arrives on a site of 
the right border of the lattice, it is removed at once. In 




FIG. 1: (Color online) 3D view of the first lattice geome- 
try we used to test our results: a cylinder (lattice with peri- 
odic boundary conditions in one direction) connected to a full 
reservoir of cells (a source) on its left end and to an empty 
reservoir of cells (a sink) on its right end. Here, the system is 
made of riL -I- 1 = 16 rings of sites (including the two rings in 
the reservoirs), and each ring is made of njj = 16 sites. After 
some transient regime, a steady current of cells establishes 
along the cylinder. A typical configuration of the cells in the 
steady state for interaction parameter p = 0.8 is shown. 



this geometry, if we let the system evolve starting from 
a random configuration, it relaxes to a steady state with 
a permanent current where gains and losses of cells from 
and to the two reservoirs compensate in average (but the 
instantaneous total number of cells still fluctuates as time 
advances). On an infinite lattice, the steady state might 
be reached only after an infinite duration. In practice, we 
simulate in parallel two families of independent systems, 
one where the site occupation probability is initially close 
to one and one where it is initially close to zero, and we 
stop our simulation when the values of the macroscopic 
observables (averaged over the systems of each family) 
are undistinguishable up to the error bars. The number 
of systems in each family is chosen so that the error bars 
are shorter than the precision that we request in advance. 

In the steady state, the concentration profile can be 
computed in the macroscopic limit thanks to Pick's 
law P?)) or to the PDE In the stationary regime 

(where no macroscopic quantity depends on time), any 
initial asymmetry has been "forgotten" by the system 
and p{j) does not depend on y because the system is 
translation-invariant in the y direction. Purthermore, 
p(r) is such that the current of cells is uniform: J, the 
density of cell current, is independent of x and y. Then 
Pick's law ini) reads 

{2l^)D{p{x)\d^p{x)^-j. (19) 
With the boundary conditions /9(0) = 1 and p{V) — 
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FIG. 2: ( Color online) Dependence on the interaction param- 
eter p of the stationary profiles p{x/L) of the site occupation 
probabiUty along the cylinder between the full reservoir (left, 
X — Q, site occupation probability p = 1) and the empty 
reservoir (right, x = L, p = Q) — see Fig. ([T}. Since the cell 
concentration c is proportional to the site occupation proba- 
bility p, the profiles of c would be the same. For each value 
of p, we plot the simulation results for the cellular automaton 
(circles) and our prediction from the analytic approximation 
Eq. (I21| (solid lines). From bottom to top, p = 0.01, 0.05, 
0.2, 0.4, 0.7 and 1. The error bars are smaller than the circles' 
size. The cylinder is made of 64 x 64 sites except for p — 0.01 
and 0.05 where there are 512 x 512 sites. 



and the expression for D{p) above, we find in 2D 

j = (l+p)/(6V3L) (20) 

and 

x{p) = {l + p[{2p-l)p{p~3)-3{l-p)]/{p+l)}L. (21) 

Then one can get the curve for p{x) by plotting the para- 
metric curve {x{p),p). 

In Fig. [2l we plotted the stationary profile of the site 
occupation probability cell p on the cylinder, between the 
two reservoirs, as obtained from simulations of the cel- 
lular automaton and from the macroscopic PDE. Notice 
the effect of nonlinear diffusion: for p = 1/2 (not shown 
in Fig. [5]), the profile is linear as predicted by the usual 
Fourier or Fick law, but, for other values of p, interac- 
tions between cells make it more complicated. For most 
values of p, there is an excellent agreement between the 
simulation results and the analytical approximation for 
all lattice sizes we tried (starting from L — 16). How- 
ever, we can distinguish two situations in which the mi- 
croscopic and the macroscopic models disagree. The first 
one is related to a boundary effect: at the vicinity of the 
empty reservoir when p = 1 (respectively at the vicinity 
of the full reservoir when p — 0), a clear discrepancy of 




FIG. 3: (Color online) Finite-size effects for the stationary 
profile of the cell concentration between the two reservoirs. 
Main curve: dots: profile p(x/L) of the site occupation 
probability from simulation results for p = 0.01 and system 
sizes ul ~ uh — 1 = 15, 31, 63, 127, 255, and 511 (from 
bottom to top). It shows a strong finite-size dependence (for 
p = 0.05, not shown, a weaker finite-size dependence is found). 
The solid line is the analytic approximation, clearly in dis- 
agreement with the simulation results. Bottom inset: for 
p=0.01, the extrapolation of the occupation probability p at 
position a; = I//4 to an infinite lattice shows that the discrep- 
ancy between the simulation data (dots) and the analytical 
approximation (dashed line) is not due to finite-size effects. 
Indeed, there is convergence toward a value of the site occu- 
pation probability that is clearly below the analytic approxi- 
mation. Top inset: profile p{x/L) for p = 1, with numerical 
data for the sizes ul = 15, 31, and 63 only (dots) and an- 
alytic prediction (solid line). The finite-size effects are very 
weak and there is a good agreement between simulation data 
for large systems and the analytic prediction. This is true 
more generally for all values of p away from (say, larger 
than 0.1). 



the average occupation probability p of the last-but-one 
lattice site (respectively of the first few lattice sites) is 
seen between the cellular automaton results and the cor- 
responding value from our analytic formula. The second 
situation is what happens when p gets close to 0; see, 
e.g., p = 0.01 and 0.05 in Fig. [2l Here, the concentra- 
tion profiles of the finite cellular automata show a strong 
finite-size dependence (main plot of Fig. [31 compared to 
the top inset of Fig. [3]), and they seem to converge to 
some shape that is clearly different from the profile pre- 
dicted from the PDE. To exclude the possibility that the 
simulation results agree with the analytic approximation 
at very large sizes, but not at the sizes we have simulated, 
we plot in the bottom inset of Fig. [3] an extrapolation 
to infinite systems: it shows that the discrepancy should 
persist even for huge systems. These phenomena are con- 
firmed in Fig. m where the permanent current of cells 
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FIG. 4: (Color online) Stationary current of cell density be- 
tween the full and the empty reservoirs, across the cylinder 
(geometry defined in Fig. [TJ, as a function of the interaction 
parameter p. To be able to compare different system sizes 
(bottom to top : nH — 1 = nz, = 15, 63 and 255), we actually 
plot the current density j multiplied by the system length 
L (dots). For ul ~ 255, only the points that differ signifi- 
cantly from the other sizes are represented. The error bars are 
smaller than the symbol's sizes, except for hl = 255, where 
they are drawn. Solid line: analytic approximation, Eq. (|20p 
times L. Like for the profile of the cell concentration, there 
is a large disagreement between the analytic approximation 
and the simulation results when p is close to and a small 
disagreement when p is close to 1. 



across the cylinder is plotted as a function of p: one sees 
a good agreement for most values of p, with a slight devi- 
ation and practically no finite-size effects around p — I, 
and some more important deviation, together with strong 
finite-size effects, around p = 0. This will be understood 
in the next subsection. 



D. Discussion of our approximations: The 
occurrence of correlations 

Our approximate computation of the hydrodynamic 
limit fails if the correlations between the occupation of 
neighboring sites are not negligible or if the continuous 
space limit has singularities (e.g., an infinite gradient) 
that cannot exist on the lattice. 



1. Singularities in continuous space 

Singularities are easy to predict and to address; they 
could happen for our model since the diffusion coefficient 
vanishes (for p = 1 and p = on one hand and for p — 



and p = 1 on the other). Indeed, Pick's law reads 

j{f,t)^-D[p{f,t)]Wp{f,t), (22) 

and at a point where D vanishes while j is finite (which 
happens for p = 1 or between the two reservoirs of 
the previous subsection), the concentration gradient is 
infinite. To get the correct, physical behavior if necessary, 
it may suffice to reintroduce two or a few lattice sites 
around the point where the singularity appears and to 
solve separately the discrete equations (obtained after 
averaging but before the Taylor expansion) on these sites, 
and the PDE on the rest of the lattice. 

In our case, this simple procedure does not enable one 
to quantitatively (or even qualitatively) understand the 
discrepancies between the automaton and the approxi- 
mate results close to the reservoirs. This is because these 
discrepancies have to do with correlations. 



2. Correlations 

If i and j are two sites, let us consider the so-called 
connected correlation function [s^ of i and j, 

{Mt)pjit))c {[Mt) - (MtMpjit) - (pjim) 

= {p.,[t)p,{t))-{p,{t)){p,{t)). (23) 

This quantity is zero if the statistical distribution of the 
occupation numbers of the sites i and j is the same as if 
the sites i and j would be filled independently at ran- 
dom with given mean occupation probabilities. More 
precisely, in the case of a finite lattice of N sites with 
a fixed total number of cells, this definition must be re- 
placed with 

{p^{t)pMc {[p^{t)~{pm)][pAiy{NpAt)~^)l{N~l)]) 

(24) 

if one wants to detect correlations due to the interactions 
between cells that could be "hidden" behind the contri- 
bution —p{\ — p)/{N — 1) due to the constraint that the 
total number of cells is fixed. 

Our approximate computation of the hydrodynamic 
limit makes use of an assumption of statistical indepen- 
dence ("all connected correlations functions are zero"), 
and we found in our simulations that the correlations are 
indeed small (they vanish as L — > -|-oo) in almost all 
cases. 

A situation in which we can compute exactly the corre- 
lations is the equilibrium state of a translation invariant 
lattice. Since, in our stochastic automaton, any move- 
ment of a cell may be reversed and since, if the lattice 
is translation invariant, any transition from one config- 
uration of the cells on the lattice to another has the 
same rate, l/(6iV), the detailed balance condition [4§| 
is fulfilled, the equilibrium state exists, and the proba- 
bility distribution of the configurations is uniform. Then 
the connected correlation functions (with the definition 
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FIG. 5: (Color online) Schematic explanation of the corre- 
lations that appear for p close to near the full reservoir 
(left) and for p close to 1 near the empty reservoir (right). 
Left: when two cells sit next to one another close to the full 
reservoir and if p is close to 0, they prevent one another from 
moving in one of the two directions that would be allowed 
if they were alone (black arrows represent permitted moves). 
Thus accompanied cells stay longer close to the reservoir than 
single ones, and positive correlations occur (see text). Right: 
when one cell sits close to the empty reservoir and if p is 
close to 1, it cannot move unless it has a neighboring cell. 
Thus single cells stay longer close to the empty reservoir than 
accompanied ones, and negative correlations occur (see text). 



above, relevant to the case of conserved total cell num- 
ber) between any two sites are zero. We simulated the 
cellular automaton in this situation and we checked that 
our numerical data are compatible with vanishing corre- 
lations (up to the error bars) for several values of p and 
of the cell concentration; this is a check of correctness 
of the simulation program. In this case, our analytical 
approximation is exact, but the result is of course trivial 
{p is uniform). 

Let us come back to the nonequilibrium steady state 
between the two reservoirs. In such a case, it is 
well known that correlations do exist, can be quanti- 
tatively iniportant, and, above all, are generically long- 
range |68l. 169(1 . For our automaton, we found numerically 
that, for p close to 1 (0) there exists a significant negative 
(positive) correlation of the occupation probability of ad- 
jacent sites in the last (first) row of sites, parallel to the 
reservoir. This is why the macroscopic limit disagrees, 
in these situations, from the simulation. In other situa- 
tions, they decrease at least as fast as 1/L as the system 
length L goes to -t-oo. Since studying these correlations 
in more detail would go beyond the scope of the present 
article, we limit ourselves here to explaining qualitatively 
what happens close to the boundaries. The correlation 
there can be interpreted using the microscopic rules (see 
Fig. [5]). If p = 1, when two cells are close to the empty, 
absorbing reservoir (Fig. O right), sooner or later one of 
the two cells will fall into the reservoir, leaving the other 
one alone. This other cell will not move until another 



cell comes on a neighboring site. Therefore, isolated cells 
stay longer close to the reservoir than accompanied cells, 
and there is more chance to find an empty site close to a 
cell than one would expect if cells were distributed at ran- 
dom, hence the negative correlations. Far away from the 
reservoir, the cell concentration is higher and there is al- 
ways some cell to "unblock" an isolated cell, so the same 
phenomenon cannot be observed. Conversely, if p = 0, 
when two cells are close to the full reservoir (Fig. [5l left) , 
they are in some sense in each other's way and they may 
progress away from the reservoir only in one of the two 
directions they could use if they were alone, the other 
direction being forbidden by the kinetic rule of the au- 
tomaton since they would stay in contact with a cell they 
were already in contact with. If one could consider only 
these cells on the nearest line of sites from the reservoir 
and forget about the cells that lie on the next-to-nearest 
line of sites, one could conclude that accompanied cells 
stay longer close to the reservoir than single cells, that 
there is more chance to find another cell close to a cell 
than one would expect if cells were distributed at ran- 
dom, and that the correlation is positive. This is not 
so simple because the next-to-nearest line of sites has 
many cells that can compensate this effect, and we do 
observe that correlations extend over a few lines of sites 
away from the full reservoir, but the qualitative idea is 
correct. Far away from the full reservoir, the site occupa- 
tion probability is not close to 1 and there is always some 
hole to "unblock" a cluster of cells, hence no such corre- 
lation is observed. We expect that taking into account 
these correlations into the analytic computation would 
improve the agremeent between simulations and analytic 
formulas. 

Now let us come to the disagreement of the simulated 
and predicted concentrations of cells in the whole system 
when p comes close to zero. Our interpretation is the 
following. We know from a preceding discussion that, 
at p = 0, the cellular automaton is a cooperative KCM 
for which Fick's law is obeyed only above some spatial 
scale that diverges with the cell concentration. Since, 
between the two reservoirs, all concentrations can be ob- 
served, Fick's law might never hold even in arbitrarily 
large systems, but this is beyond the scope of the present 
paper. If p is close to zero but strictly positive, we expect 
that the threshold spatial scale for effectively Brownian 
diffusion should be large (diverging with p) but finite at 
all values of the cell concentration. Therefore, and con- 
trary to what happens for larger values of p where the 
threshold scale must be comparable to the lattice spacing, 
there is a whole range of system sizes where the Brownian 
regime is not yet reached and the concentration profile 
depends strongly on the system size. Above this scale, 
diffusion sets in but, as a consequence of the cooperative 
nature of the kinetics, correlations always play a signif- 
icant role and our evaluation of the diffusion constant, 
which neglects them, is irrelevant, hence the disagree- 
ment observed even at large system sizes between the 
predicted profile and the actual one. 
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IV. RESULTS FOR CELLS DIFFUSING OUT OF 
A SPHEROID 

In this section, we come back to the radial geometry, 
where cells diffuse away from a center (cf. Sec. II B), in 
order to compare the results from the macroscopic model 
to experimental results. We use the expression of the dif- 
fusion coefficient obtained in Sec. Ill B 2, cf. Eq. (|14p . for 
a 2D system on a hexagonal lattice. We calculate numer- 
ically the solution of the nonlinear diffusion equation in 
a radial geometry by using a fully implicit method where 
the spatial derivatives in the diffusion equation are evalu- 
ated at time step t+1 (for the sake of stability) [zO]. The 
site occupation probability p is constrained to be equal 
to 1 in the center at all times. 

A comparison between the density profiles from the au- 
tomaton and from the numerical solution of the diffusion 
equation at the same time is presented in Fig.[Bl As in the 
case of diffusion between a full and an empty reservoir, 
the agreement is very good except for p — 0.05, close to 
the full reservoir, where the numerical solution does not 
agree with the density profile of the automaton as nicely 
as in the cases of larger p. This effect is due to corre- 
lations between close neighbors which are not negligible 
when p is close to (cf. Sec. Ill D 2). In general, the 
agreement between both profiles is nice. This is a supple- 
mentary check of our derivation of a nonlinear equation 
from the automaton. 
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FIG. 6; Density profiles from tlie automaton (dotted line) 
and from the numerical solution of the diffusion equation (full 
line) for p = 1, 0.5, and 0.05. Error bars for the automaton 
density profile are smaller than the markers and thus are not 
represented. 

This good agreement allows us to make a direct com- 
parison between the numerical solution of the nonlinear 
diffusion equation and experimental cell densities. The 
experiments consisted in placing a spheroid of glioma 
cells (GL15) on a collagen substrate and following the 



evolution of the migration pattern during a few days. 
A detailed protocol has been described in [s^]. Exper- 
iments have been performed by Christov at the Henri 
Mondor Hospital. Figure [7] is a photograph of an ex- 
perimental spheroid, 36 h after cells started migrating. 
In order to compare the experimental density profiles at 
different times of evolution (cf. the full lines in Fig. [5]) 
and the numerical solution of the diffusion equation, we 
rewrite Eq. (fT4|) as 

D{p) = D., [2(1 -p) + 4(2p - l)p - 2{2p - l)p\ (25) 

where is the value of D when P ~ \ (i-e., with- 
out interactions between cells). Note that = i when 
the time step of the numerical solution coincides with 
that of the automaton (cf. Sec. Ill B 2). For compar- 
ison with experimental data, the numerical time step is 
chosen equal to 1 h (and thus I?* can be different from 
i). The spatial length step is set to 35.2 /im as in the 
automaton, which corresponds to a characteristic size of 
cells. The diffusion equation was integrated for p — 0.95. 
The choice p — I induces sharp edges at low densities 
that do not correspond to our experimental profiles. 

We find that, for the coefficient = 1240 ± 

100 /um^h~^, numerical solutions at 12 and 36 h agree 
well with the experimental data at the same times, cf. 
Fig. [5] One can notice that this value falls within the 
range of values of the diffusion coefficient values ob- 
tained in vivo '71j (between 4 and 88 mm^/d = 450 — 
10000 pni^ h^^), although this comparison must be taken 
with a grain of salt due to the different nature of the sys- 
tems. 




FIG. 7: Experimental pattern of migration after 36 h. The 
real size of the image is 1.90 x 1.77 mm. 
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FIG. 8: Comparison between the numerical solution of the 
nonlinear diffusion equation (full lines) and experimental 
data, at t = 12 h (crosses) and t = 36 h (circles). D* — 
1240 ± 100 fim'^ h-^ with a = 35.2 ^m. 

For this type of in vitro experiments, where cells are 
followed individually, the automaton remains the best 
way to model the system: time is counted as the number 
of cells that have left the spheroid, the lattice step corre- 
sponds to the mean size of a cell, and interaction between 
cells is represented by rules defined between neighboring 
cells. The only parameter to be fixed is p. On the con- 
trary, for real tumors where the number of cells is of the 
order of 10^, the automaton becomes unmanageable, and 
in this case, the approach through a diffusion equation is 
more efhcient. 



V. CONCLUSION AND OUTLOOK 

In this paper we set out to establish a link between 
a microscopic description of tumor cell migration with 
contact interactions (which we developed in 35]) and a 
macroscopic ones, based on diffusion-type equations used 
in most phcnomenological models of tumor growth. 

A microscopic description has obvious advantages. 
First, it allows a straightforward visualisation of the re- 
sults: one "sees" the cells migrating. Moreover it can be 
directly compared to experiments and thus reveal special 
features of the migration like, for instance, cell attrac- 
tion which we postulated in (35| in order to interpret the 
experimental cell distribution. It allows predictions for 
the trajectories of individual cells that can be observed 
in time-lapse microscopy experiments [3^ [72l . [73| . How- 
ever the microscopic approach is not without its draw- 
backs. In all implementations we presented we have lim- 
ited ourselves to a two-dimensional geometry. While a 
three-dimensional extension is, in principle, feasible it 
would lead to a substantial complication of the model. 
Moreover, even in a 2-D setting, it is not computation- 
ally feasible to accomodate the millions of cells present 



in a tumor. Our typical migration simulations involve at 
best a few thousand cells. 

A macroscopic treatment on the other hand deals with 
mean quantities, like densities. Thus limitations such as 
the number of elementary entities (cells, in the present 
work) simply disappear. However the trace of micro- 
scopic interactions equally disappears unless it is explic- 
itly coded into the macroscopic equations. This is a very 
delicate matter: the present work is dedicated to bridging 
just this gap between the microscopic and macroscopic 
approaches. 

Our starting point was the microscopic description of 
cell migration which was validated through a direct com- 
parison to experiment. Starting from a cellular automa- 
ton (with a given geometry and update rules) we have 
derived a nonlinear diffusion equation by taking the ad- 
equate continuous limits. Our main finding was that the 
interactions between cells modify the diffusion process 
inducing nonlinearities. 

Such nonlinearities have been observed in models of 
ecosystems where individuals tend to migrate faster in 
overcrowded regions [t^, IzE]- This mechanism (fleeing 
overcrowded regions) has also been the subject of spec- 
ulation in the context of cell migration [76^, but, in our 
experimental situation, it is irrelevant: on one hand, cells 
in overcrowded regions cannot migrate fast, even if they 
want to, because they are fully packed (thus blocked); 
on the other hand it has been shown experimentally that 
the nonlinearity is due to attractive interactions between 
cells (when interactions are suppressed, the nonlinearity 
diminishes) [Ab]. 

While the present model is most interesting, it is fair to 
point out the precise assumptions that entered its deriva- 
tion. We were based on a specific cellular automaton: 
many more do exist, the limit being only one's imagina- 
tion. Of course for the problem at hand it was essential 
to have a geometry allowing every site to possess a suffi- 
ciently large number of neighbours so as not to artificially 
hinder the cell motion. Moreover the update rules were 
chosen in fine through a comparison to experiment. But 
it would be interesting to see how robust the nonlinear 
diffusion coefficient is with respect to changes of the ge- 
ometry: introducing a random deformation of the lattice 
or mimicking the migration of cancer cells on a susbtrate 
of astocytes in the spirit of [45]. Moreover, if one for- 
gets about the application to tumor cell migration, there 
may exist several interesting cases of cellular automaton 
geometry and kinetic rules leading to nonlinear diffusion 
equations worth studying, in addition to the kinetically 
constrained models we have seen and that have their own 
interest for the physics of glasses. We intend to come 
back to these questions in some future publications. 

Finally, an incontrovertible generalization of the 
present work would be that of three-dimensional mod- 
els. While the extension of microscopic models to three 
dimensions may be computationally overwhelming, the 
solution of three-dimensional diffusion equations, be they 
nonlinear, does not present particular difficulties. Thus 
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once the automaton rules (mimicking the cell-cell inter- Acknowledgments 
action) have been used for the derivation of a diffusion 
equation like Eq. (|13p with coefficient (fT5| . we can pro- 
ceed to the study of the macroscopic model, coupling dif- 
fusion to proliferation (plus, perhaps, other macroscop- 
ically described effects) for the modeling of real-life tu- 
mors (in particular, glioblastomas). This is a path we We acknowledge financial support from the "Comite 
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